altman-demetrio-tarasenko Stockholm 2023 replication

A structured comment on Cavaillé and Ferwerda (2022)

Micah Altman¹² (MIT) escience@mit.edu, Hans Gaebler (Harvard) jgaebler@fas.harvard.edu, Georgy Tarasenko (Cornell) gt298@cornell.edu

Abstract

Cavaillé and Ferwerda (Cavaillé and Ferwerda 2023) examine the effect of granting immigrants access to the welfare state in Austria during 2006 on the support for far-right parties. In their preferred specification, they find that – the interaction of the proportion of non-eu residents in a municipality with the proportion of population in that municipality has an effect on the change in vote share for far right parties of +0.67 (p < .05) supporting the main substantive claim that increasing welfare benefits to immigrants increases support for far-right parties. Moderating this, based on a subset of data from Vienna, they find that – percentage of population in rental housing and the percentage of population in public housing are both substantive and significant predictors of change in far-right party vote share with magnitude of .03 and .09 (respectively) and p-values of < .05. Based on this estmate they conclude that support was increased by the the high proportion of public housing beneficiaries or low-end rentals within districts.

As part of a collaborative iniatitve, we conducted a time-bounded replication study, following a standardized template. We first reproduce the paper’s main and secondary analyses, and subject these to computational robustness tests.

We then conduct several types of conceptual replication – examining the robustness of results to selection of covariates, reweighting, and the sensitivity of results to resampling and outliers.

We find that the research is computationally replicable but the evidence for the main causal claims are potentially overstated as it sensitive to a small set of observations, does not explain most of the variation in outcomes, and does not clearly out-compete alternative specifications.

Introduction

Cavaillé and Ferwerda (Cavaillé and Ferwerda 2023) tested the effect of granting immigrants access to the welfare state in Austria during 2006 on the support for far-right parties. We focus on their two main claims, as described in the abstract (pg 19):

we show that the reform increased support for far-right parties with welfare chauvinist platforms. Electoral ward data suggest that this response was concentrated in districts with a high proportion of public housing beneficiaries or low-end rentals. Our findings provide novel evidence that distributional conflict can accelerate the rise of far-right parties in countries with substantial in-kind welfare programs

In their preferred specification, using OLS with clustered standard errors, they find that – the interaction of the proportion of non-eu residents in a municipality with the proportion of population in that municipality has an effect on the change in vote share for far right parties of +0.67 (p < .05) (pg 26, Table 1, col 1) supporting the main substantive claim that increasing welfare benefits to immigrants increases support for far-right parties (pg. 19). Moderating this, based on a subset of data from Vienna, they find that – percentage of population in rental housing and the percentage of population in public housing are both substantive and significant predictors of change in far-right party vote share with magnitude of .03 and .09 (respectively) and p-values of < .05. (pg 30, Table 2, col 1). Based on this estimate they conclude that support was accelerated by the the high proportion of public housing beneficiaries or low-end rentals within districts. (pg 19.)

Analysis Protocol and Context

We created this replication report as part of an intensive single-day replication effort as part of the ‘Replication Games’ a multi-analyst research initiative of the Institute for Replication. We pre-selected the study from a short list of candidates provided by the organizers, and applied the provided template and nomenclature.

The protocol for the overall initiative was pre-registered by the principle investigators of the effort. ³ The protocol does not specify the details of interaction with or instruction to the analysts, which were later provided to us as slides and email. In summary, analysts were asked to: First, approximately 1-2 weeks prior to the games familiarize ourselves with the article to be replicated, and develop an overall plan for replication day. 2. During replication day, conduct a 1-day (8-hour) intensive replication of the article. (Analysts were instructed that it was permissible to continue after this period, but encouraged to focus their effort and planning on this day.) 3. Over approximately the next 6-8 weeks write up the results of the replication-day exercise using a standard template.

Our team’s work was consistent with these instructions. Approximately one day of effort each was allocated by the team to planning, replicating, and writing ( for a total of 9 FTE full-time-equivalent days of effort).

Subsequently, we were asked to review a response from the author to our report. We allocated approximately two days to this activity – which resulted in some additions to the text of this report to clarify its interpretation.

Overview of Analysis

We obtained the published replication data set (Cavaille 2023) from the Harvard Dataverse archive, where it had been deposited for public use by the authors. The replication data set included the final processed data used to produce the paper results, along with R code to replicate figures and tables. The data was sufficient to check computational reproducibility but posed challenges for conceptual reproducibility because it provided neither copies of the source data, nor links to or citations of the original data sources (which were described generally). Further although the replication data did include additional measures not used in publication, these were largely undocumented, and hence difficult to reliably interpret.

We first reproduced the paper’s main and secondary analyses, and subject these to computational robustness tests. The computational reproductions showed both the main and secondary results to be reproducible. As an unintended consequence reproducing the results revealed the substantively poor fit of these models: the primary model has an R-squared of .0356. This finding informed further replication analysis.

We then conducted several types of conceptual replication – examining the robustness of results to selection of covariates, reweighting, and the sensitivity of results to resampling and outliers.

Prompted by the unexpectedly poor fit of the primary model, we performed a conceptual replication that applied the authors’ linear model with clustered errors with alternate covariates, and to models that did not include the interaction effects that are central to the papers proposed causal explanation. This reveals that the addition of the model interaction term does not a provide substantial improvement in predictive power over plausible alternative models.

The authors suggest, in their first analysis of nation-wide election data, that voters in municipalities with high proportions of residents living in public housing as well as comparatively high proportions of third-country nationals voted more for far-right parties as a result of the demand shock on public housing. We undertake a conceptual replication of this result using their data on the city of Vienna. In particular, we look at census tracts—as a proxy for neighborhood—to see if tracts with a higher proportion of third-country nationals and higher proportions of individuals living in public housing tend to vote for far-right parties at higher rates. We find that implementing this conceptual replication increases the magnitude of the main point estimate for the interaction between these two factors, but the estimate is no longer statistically significant at the 5% level.

The authors fit the models in both of their main analyses at the level of administrative units—i.e., all municipalities are equally weighted in their nation-wide analysis, and all tracts are equally weighted in their analysis of Vienna. However, it is unclear if this weighting is substantively appropriate. To test the robustness of their conclusions to this, we reweight their model so that each administrative unit has weight equal to the number of voters it contains. Implementing this robustness check in their primary analysis increases the magnitude of their main point estimate for the the interaction of the proportion of residents living in public housing with the proportion of third-country nationals residing in a municipality, which remains significant at the 5% level. Implementing this robustness check in their secondary analysis has no effect on the magnitude or the statistical significance of the main point estimates.

Next, to better understand the authors’ main claims, how well their models fit the data, and to search for potential anomalies not evident in low-dimensional summaries, we attempted to visually replicate their primary and secondary analyses. Based on this graphical exploration, we carried out an outlier analysis.

We assessed the robustness of the findings from the Austrian sample by examining potential outliers in the key variables. We observed that the previously significant interaction effect between the share of non-EU residents and the share of public housing loses its significance when we exclude abnormally large values of key variables. Covariate valance analysis highlights substantial differences between observations with outliers and those without. When we reevaluate the main argument using the sample of outliers, we find that all the initial effects become notably stronger in both magnitude and significance. This may suggest that the authors’ primary argument may not generalize to all Austrian districts within the sample.

Finally, given the substantive importance placed on the interaction term in the authors’ primary analysis, as well as the generally low goodness-of-fit of the models considered, we evaluated the main models of the primary and secondary analyses using out-of-sample tests of predictive accuracy. We find in both cases that the models with and without the interaction term have essentially the same mean squared error, suggesting that the substantive interpretation of the interaction effect may not be appropriate.

Reproducibility

Simple Direct Reproduction (Calibration)

As a baseline check of the model, and our understanding of it we conducted a simple computational replication, using the authors’ supplied data and code, and compared these to published results.

Reproduction of Table 1 - Model 1 (Primary Result)

Loading required namespace: haven

Loading required namespace: estimatr

Reproduction of Table 2 - Model 1 (Secondary Result)

vienna_authors.df <- 
  haven::read_dta("authors replication materials/vienna_final.dta")

m2.formula<- formula("dv ~  (pctrental + pctpublic_w_zsp)")

results_m2_computational.lmr <-
  estimatr::lm_robust(m2.formula,   data = vienna_authors.df, 
                      clusters = tract_key)

Warning in eval(quote({: Some observations have missingness in the cluster
variable(s) but not in the outcome or covariates. These observations have been
dropped.

tidy_m2_original.df <-
  structure(
    list(
      term = c(
        "(Intercept)",
        "pctrental",
        "pctpublic_w_zsp"
      ),
      estimate = c(0.04, 0.03, 0.09),
      std.error = c(0.01, 0.01, 0.01),
      p.value = c(0.05, 0.05,  0.05),
      repl_id = c("original", "original", "original")
    ),
    row.names = c(NA,-3L),
    class = "data.frame"
  )

repro_m2.df <- tidy_results(list(comp=results_m2_computational.lmr)) %>% 
   bind_rows(tidy_m2_original.df)

repro_m2_summary.df <- tidy_summary(list(comp=results_m2_computational.lmr))

repro_m2.df %>% group_by(repl_id) %>% 
      fmtTB(caption="coefficients")

repro_m2_summary.df %>% group_by(repl_id) %>%
      fmtTB(caption="summaries")

Summary - Direct Reproducibility

The computational reproductions showed both the main and secondary results to be reproducible.

Computational robustness

Generally, even with a fixed statistical model family and specification results may vary with the estimation algorithm used and specific software’s implementation of it. (Altman, Gill, and McDonald 2004) We evaluate the computational robustness of the model by using alternative algorithms and implementations.

Primary Result

requireNamespace("arm")

Loading required namespace: arm

results_m1_statlm.lm <-
  stats::lm(m1.formula, data = austria_authors.df)
results_m1_statglm.glm <-
  stats::glm(m1.formula, data = austria_authors.df)
results_m1_statbayeslm.glm <-
  arm::bayesglm(m1.formula, data = austria_authors.df, 
                prior.scale=Inf, prior.df=Inf)

#Note: could also add nls, bayesglm, mle2 -- would require re-expressing current formula in different syntax, and renaming results matrix

ml.ls <- list(lmrobust = results_m1_computational.lmr,
              lm=results_m1_statlm.lm,
              glm=results_m1_statglm.glm,
              bayes=results_m1_statbayeslm.glm)

alt_est_m1.df <- tidy_results(ml.ls)

Warning: The `tidy()` method for objects of class `bayesglm` is not maintained by the broom team, and is only supported through the `glm` tidier method. Please be cautious in interpreting and reporting broom output.

This warning is displayed once per session.

alt_est_m1_summary.df <- tidy_summary(ml.ls)

alt_est_m1.df  %>% group_by(repl_id) %>% 
        fmtTB(caption="coefficients")

alt_est_m1_summary.df %>% group_by(repl_id) %>% 
        fmtTB(caption="summaries")

rm(ml.ls,results_m1_statlm.lm,results_m1_statglm.glm,results_m1_statbayeslm.glm)

Secondary Result

requireNamespace("arm")

results_m2_statlm.lm <-
  stats::lm(m2.formula, data = vienna_authors.df)
results_m2_statglm.glm <-
  stats::glm(m2.formula, data = vienna_authors.df)
results_m2_statbayeslm.glm <-
  arm::bayesglm(m2.formula, data = vienna_authors.df, 
                prior.scale=Inf, prior.df=Inf)

#Note: could also add nls, bayesglm, mle2 -- would require re-expressing current formula in different syntax, and renaming results matrix

ml.ls <- list(lmrobust = results_m2_computational.lmr,
              lm=results_m2_statlm.lm,
              glm=results_m2_statglm.glm,
              bayes=results_m2_statbayeslm.glm)

alt_est_m2.df <- tidy_results(ml.ls)
alt_est_m2_summary.df <- tidy_summary(ml.ls)

alt_est_m2.df  %>% group_by(repl_id) %>% 
        fmtTB(caption="coefficients")

alt_est_m2_summary.df %>% group_by(repl_id) %>%
        fmtTB(caption="summaries")

rm(ml.ls,results_m2_statlm.lm,results_m2_statglm.glm,results_m2_statbayeslm.glm)

Summary - Computation Robustness of Reproducibility

The computational reproductions showed both the main and secondary results to be robust to alternative choices of algorithm and software implementation.

Replication:

Conceptual Replication: Alternate OLS Specification (covariate robustness)

Presenting multiple models - as opposed to pure null-hypothesis-significance-tests of a main result - is common in political science, where much analysis is observational and nature, and there are often plausible competing causal explanations. Consistent with this, the original article presents results for multiple alternative models of the theorized interaction effect. In particular, Table 1 of the original paper presents eight different variants of the interactive model, with associated coefficient estimates. However, neither the paper nor the supplementary materials provided goodness-of-fit criteria for the models; nor explicitly compared the models to competing models in which interactions were absent.

As a routine part of replicating the models, a number of standard diagnostic measures of fit were produced. Examing these revealed, unexpectedly, the the primary model had an R-squared of .0356 – which is presumptively indicative of a poor model fit, suggests that the model is mispecified. This finding informed the replication analysis.

Prompted by the unexpectedly poor fit of the primary model, we performed a conceptual replication that applied the authors’ linear model with clustered errors with alternate covariates, and to models that did not include the interaction effects that are central to the papers proposed causal explanation.

Alternate Covariates & Model 1

#NOTE: could also explore  alternate methods for computing robust standard errors (e.g. sensemaker, lmtest, sandwich)

results_m1_clustered.lmr <-  estimatr::lm_robust(m1.formula, data = austria_authors.df, clusters = bezirk)

results_m1_interactionsonly.lmr <-
  estimatr::lm_robust(formula("d_rr_06 ~ dv_pop_01:pct_noneu_06 -dv_pop_01 -pct_noneu_06"),
                      data = austria_authors.df)

results_m1_kitchensink.lmr <-
  estimatr::lm_robust(
    formula(
      "d_rr_06 ~ dv_pop_01:pct_noneu_06 +educ_tertiary +avg_income +lab_pct_manufact_01 +lab_pct_unemp +welfare_cap_06 +health_cap_06 +education_cap_06 +foreignborn_delta+citizen_eu_growth_pct +vacancy_01_public"
    ),
    data = austria_authors.df
  )

results_m1_kitchensinkmain.lmr <-
  estimatr::lm_robust(
    formula(
      "d_rr_06 ~ dv_pop_01+pct_noneu_06 +educ_tertiary +avg_income +lab_pct_manufact_01 +lab_pct_unemp +welfare_cap_06 +health_cap_06 +education_cap_06 +foreignborn_delta+citizen_eu_growth_pct +vacancy_01_public"
    ),
    data = austria_authors.df
  )

# vacancy has high missingbess
results_m1_kitchensink_novacancy.lmr <-
  estimatr::lm_robust(
    formula(
      "d_rr_06 ~ dv_pop_01:pct_noneu_06 +educ_tertiary +avg_income +lab_pct_manufact_01 +lab_pct_unemp +welfare_cap_06 +health_cap_06 +education_cap_06 +foreignborn_delta+citizen_eu_growth_pct "
    ),
    data = austria_authors.df
  )

results_m1_kitchensinkmain_novacancy.lmr <-
  estimatr::lm_robust(
    formula(
      "d_rr_06 ~ dv_pop_01+pct_noneu_06 +educ_tertiary +avg_income +lab_pct_manufact_01 +lab_pct_unemp +welfare_cap_06 +health_cap_06 +education_cap_06 +foreignborn_delta+citizen_eu_growth_pct "
    ),
    data = austria_authors.df
  )

results_m1_mainonly.lmr <-
  estimatr::lm_robust(formula("d_rr_06 ~ dv_pop_01 + pct_noneu_06 "), data = austria_authors.df)

results_m1_euonly.lmr <-
  estimatr::lm_robust(formula("d_rr_06 ~ pct_noneu_06 "), data = austria_authors.df)

results_m1_poponly.lmr <-
  estimatr::lm_robust(formula("d_rr_06 ~ dv_pop_01 "), data = austria_authors.df)

results_m1_null.lmr <-
  estimatr::lm_robust(formula("d_rr_06 ~ 1 "), data = austria_authors.df)

ml.ls <- list(
  author_model1 = results_m1_computational.lmr,
  author_clusterederr = results_m1_clustered.lmr,
  loaded_model = results_m1_kitchensink.lmr,
  loaded_main = results_m1_kitchensinkmain.lmr,
  loaded_nv = results_m1_kitchensink_novacancy.lmr,
  loaded_main_nv = results_m1_kitchensinkmain_novacancy.lmr,
  interaction = results_m1_interactionsonly.lmr,
  maineffects = results_m1_mainonly.lmr,
  single_eu = results_m1_euonly.lmr,
  single_pop = results_m1_poponly.lmr,
  null_model= results_m1_null.lmr
)

alt_var_m1.df <- tidy_results(ml.ls)
alt_var_m1_summary.df <- tidy_summary(ml.ls)

alt_var_m1.df  %>% group_by(repl_id) %>% 
        fmtTB(caption="coefficients")

alt_var_m1_summary.df %>% group_by(repl_id) %>% 
        fmtTB(caption="summaries")

rm(
  ml.ls,
  results_m1_clustered.lmr,
  results_m1_kitchensink.lmr,
  results_m1_interactionsonly.lmr,
  results_m1_kitchensink_novacancy.lmr,
  results_m1_kitchensinkmain_novacancy.lmr,
  results_m1_mainonly.lmr,
  results_m1_euonly.lmr,
  results_m1_poponly.lmr,
  results_m1_kitchensinkmain.lmr
)

Observe that none of the alternative models fit the data well. Moreover, models without interaction effects provide only slightly lower explanatory power, as measured by adjusted R-squared; while alternative specifications with additional variables improve predictive power.

As an addenda, after reflection on the authors’ response to the first version of this comment, we emphasize that these measures are not included because they are of substantive theoretical interest as estimands in their own right, but as the default (produced by the model) measure of fit used to compare competing explanatory models. And, of course, there are many other systematic approaches to model selection and model uncertainty that could be applied – AIC, BIC, Bayes Factor, evidence , etc. (Kadane and Lazar 2004)

(We include for illustration, below, some alternate measures for model comparison. Some of these –particularly the evidence ratio– are more favorable to the authors’ preferred specification. However, neither the authors’ originally published results, nor our replication compel rejection of plausible competing models. )

In this comment, we do not argue for the adoption a specific alternative model, nor for a single measure or method of model comparison. We do observe that the addition of model fit criteria, and a systematic comparison to alternative models reflecting from competing (e.g. simpler) explanation should change a reader’s assessment of the weight of the evidence in this case, and would be a useful practice for general adoption in observational social-science analyses.

#recomput using lm -- since estimatr::  does not have an AIC method

results_m1_author.lm <-  lm(m1.formula, data = austria_authors.df)

results_m1_mainonly.lm <-
  lm(formula("d_rr_06 ~ dv_pop_01 + pct_noneu_06 "), data = austria_authors.df)

results_m1_kitchensink.lm <-
  lm(
    formula(
      "d_rr_06 ~ dv_pop_01:pct_noneu_06 +educ_tertiary +avg_income +lab_pct_manufact_01 +lab_pct_unemp +welfare_cap_06 +health_cap_06 +education_cap_06 +foreignborn_delta+citizen_eu_growth_pct +vacancy_01_public"
    ),
    data = austria_authors.df
  )

results_m1_kitchensinkmain.lm <-
  lm(
    formula(
      "d_rr_06 ~ dv_pop_01+pct_noneu_06 +educ_tertiary +avg_income +lab_pct_manufact_01 +lab_pct_unemp +welfare_cap_06 +health_cap_06 +education_cap_06 +foreignborn_delta+citizen_eu_growth_pct +vacancy_01_public"
    ),
    data = austria_authors.df
  )
results_m1_kitchensink_novacancy.lm <-
  lm(
    formula(
      "d_rr_06 ~ dv_pop_01:pct_noneu_06 +educ_tertiary +avg_income +lab_pct_manufact_01 +lab_pct_unemp +welfare_cap_06 +health_cap_06 +education_cap_06 +foreignborn_delta+citizen_eu_growth_pct "
    ),
    data = austria_authors.df
  )

results_m1_kitchensinkmain_novacancy.lm <-
  lm(
    formula(
      "d_rr_06 ~ dv_pop_01+pct_noneu_06 +educ_tertiary +avg_income +lab_pct_manufact_01 +lab_pct_unemp +welfare_cap_06 +health_cap_06 +education_cap_06 +foreignborn_delta+citizen_eu_growth_pct "
    ),
    data = austria_authors.df
  )

results_m1_null.lm <-
 lm(formula("d_rr_06 ~ 1 "), data = austria_authors.df)

ml.ls <- list(
  author_model1 = results_m1_author.lm,
  loaded = results_m1_kitchensink.lm,
  loaded_main = results_m1_kitchensinkmain.lm,
  maineffects = results_m1_mainonly.lm,
  loaded_nv = results_m1_kitchensink_novacancy.lm,
  loaded_main_nv = results_m1_kitchensinkmain_novacancy.lm,
  null_model= results_m1_null.lm
)

bictab.tb <- AICcmodavg::bictab(ml.ls)
bictab.tb 

Model selection based on BIC:

                K       BIC Delta_BIC BICWt Cum.Wt      LL
author_model1   5 -10146.28      0.00  0.94   0.94 5092.57
maineffects     4 -10140.77      5.51  0.06   1.00 5085.93
null_model      2 -10098.80     47.48  0.00   1.00 5057.18
loaded_nv      12  -9674.52    471.76  0.00   1.00 4883.56
loaded_main_nv 13  -9665.73    480.55  0.00   1.00 4883.02
loaded         13  -5643.43   4502.85  0.00   1.00 2868.40
loaded_main    14  -5633.61   4512.67  0.00   1.00 2867.08

AICcmodavg::evidence(bictab.tb, model.high="author_model1", model.low="maineffects")

Evidence ratio between models 'author_model1' and 'maineffects':
15.75 

rm(
   ml.ls,
   results_m1_author.lm,
   results_m1_kitchensink.lm,
   results_m1_kitchensinkmain.lm,
    results_m1_mainonly.lm,
   results_m1_kitchensink_novacancy.lm,
    results_m1_kitchensinkmain_novacancy.lm,
   null_model= results_m1_null.lm,
   bictab.tb
 )

Alternate Covariates & Model 2

#m2.formula<- formula("dv ~  (pctrental + pctpublic_w_zsp)")
#results_m2_computational.lmr <-
#  estimatr::lm_robust(m2.formula,   data = vienna_authors.df, 
#                      clusters = tract_key)


# m1.formula  <- formula("d_rr_06 ~ dv_pop_01*pct_noneu_06")

results_m2_interaction.lmr <-
  estimatr::lm_robust(formula("dv ~ pctrental*pctpublic_w_zsp"),
                      data = vienna_authors.df,
                      clusters = tract_key)

Warning in eval(quote({: Some observations have missingness in the cluster
variable(s) but not in the outcome or covariates. These observations have been
dropped.

results_m2_m1model.lmr <-
  estimatr::lm_robust(formula("dv ~ pctpublic_w_zsp*pctforeign"),
                      data = vienna_authors.df,
                      clusters = tract_key)

results_m2_rental.lmr <-
  estimatr::lm_robust(formula("dv ~ pctrental"),
                      data = vienna_authors.df,
                      clusters = tract_key)

Warning in eval(quote({: Some observations have missingness in the cluster
variable(s) but not in the outcome or covariates. These observations have been
dropped.

results_m2_housing.lmr <-
  estimatr::lm_robust(formula("dv ~ pctpublic_w_zsp"),
                      data = vienna_authors.df,
                      clusters = tract_key)

Warning in eval(quote({: Some observations have missingness in the cluster
variable(s) but not in the outcome or covariates. These observations have been
dropped.

ml.ls <- list(author_model2 = results_m2_computational.lmr,
              interaction = results_m2_interaction.lmr,
              m1spec = results_m2_m1model.lmr,
              housing_only = results_m2_housing.lmr,
              rental_only = results_m2_rental.lmr
              )

alt_var_m2.df <- tidy_results(ml.ls)
alt_var_m2_summary.df <- tidy_summary(ml.ls)

alt_var_m2.df  %>% group_by(repl_id) %>% 
        fmtTB(caption="coefficients")

alt_var_m2_summary.df %>% group_by(repl_id) %>% 
        fmtTB(caption="summaries")

rm(
  ml.ls,
  results_m2_interaction.lmr,
  results_m2_m1model.lmr,
  results_m2_housing.lmr,
  results_m2_rental.lmr
)

Observe as with the primary finding that none of the models fit the data well by the adjusted r-square measure, and that models without interaction effects provide only slightly lower explanatory power. Further, when the original model applied to the country as a whole (as above) for the main finding is fit to this subset, the interaction term is not significant.

Proportion of third-country nationals

In light of the eye tests below, we investigate the counterintuitive relationship between the proportion of third-country nationals living in a particular ward and the change in far-right vote share. In particular, if the mechanism driving change support for far-right parties is direct competition between third-country and Austrian nationals for public housing, then, assuming individuals have a preference for continuing to live near their current residence, we would expect that voters living in areas with more third-country nationals would tend to be more supportive of far-right parties.

results_m2_pctforeign <-
  estimatr::lm_robust(dv ~ pctforeign * pctpublic_w_zsp, data = vienna_authors.df,
                      clusters = tract_key)
tidy_results(list("pctforeign" = results_m2_pctforeign)) %>%
  gt()

repl_id	term	estimate	std.error	statistic	p.value	conf.low	conf.high	df	outcome
pctforeign	(Intercept)	0.08447235	0.005546008	15.231199	1.276299e-28	0.07347881	0.09546589	107.66008	dv
pctforeign	pctforeign	-0.10894320	0.028262648	-3.854671	2.601661e-04	-0.16534381	-0.05254260	67.77636	dv
pctforeign	pctpublic_w_zsp	0.04470112	0.007389251	6.049479	1.352690e-07	0.02989163	0.05951061	54.81840	dv
pctforeign	pctforeign:pctpublic_w_zsp	0.08442037	0.052470547	1.608910	1.129095e-01	-0.02054613	0.18938687	59.72958	dv

tidy_summary(list("pctforeign" = results_m2_pctforeign)) %>%
  gt()

repl_id	r.squared	adj.r.squared	df.residual	res_var	nobs
pctforeign	0.1765565	0.1751671	1778	0.001864403	1782

results_m2_pctforeign_w_controls <-
  estimatr::lm_robust(dv ~ pctforeign * pctpublic_w_zsp + lab_pct_pensioners + educ_tertiary,
                      data = vienna_authors.df, clusters = tract_key)
tidy_results(list("pctforeign_controls" = results_m2_pctforeign_w_controls)) %>%
        fmtTB(caption="coefficients")

tidy_summary(list("pctforeign_controls" = results_m2_pctforeign_w_controls)) %>%
        fmtTB(caption="summaries")

It’s possible that these results are driven by homophily, i.e., voters residing in tracts with higher percentages of third-country nationals have political preferences that are more friendly to third-country nationals. To test this, we see if voters living in areas with few third-country nationals were also more likely to vote for far-right parties in 2002. To ensure that the controls make sense, we reverse-engineer pctforeign02.

placebo_df <- vienna_authors.df %>%
  mutate(pctforeign02 = pctforeign / (1 + pctforeign_delta))
results_m2_pctforeign_placebo <-
  estimatr::lm_robust(farright_share2002 ~ pctforeign02 * pctpublic_w_zsp, data = placebo_df,
                      clusters = tract_key)
tidy_results(list("pctforeign_placebo" = results_m2_pctforeign_placebo)) %>%
  gt()

repl_id	term	estimate	std.error	statistic	p.value	conf.low	conf.high	df	outcome
pctforeign_placebo	(Intercept)	0.0746101488	0.001937988	38.49876040	2.582460e-64	0.070768098	0.07845220	106.47809	farright_share2002
pctforeign_placebo	pctforeign02	0.0132415358	0.010541431	1.25614215	2.133133e-01	-0.007789136	0.03427221	68.79969	farright_share2002
pctforeign_placebo	pctpublic_w_zsp	0.0208641718	0.003509215	5.94553772	2.068608e-07	0.013829025	0.02789932	54.13661	farright_share2002
pctforeign_placebo	pctforeign02:pctpublic_w_zsp	0.0009347354	0.025369795	0.03684442	9.707320e-01	-0.049817589	0.05168706	59.70385	farright_share2002

tidy_summary(list("pctforeign_placebo" = results_m2_pctforeign_placebo)) %>%
  gt()

repl_id	r.squared	adj.r.squared	df.residual	res_var	nobs
pctforeign_placebo	0.07523995	0.07367962	1778	0.0005060312	1782

results_m2_pctforeign_w_controls_placebo <-
  estimatr::lm_robust(farright_share2002 ~ pctforeign02 * pctpublic_w_zsp + lab_pct_pensioners + educ_tertiary,
                      data = placebo_df, clusters = tract_key)
tidy_results(list("pctforeign_controls_placebo" = results_m2_pctforeign_w_controls_placebo)) %>%
        fmtTB(caption="coefficients")

tidy_summary(list("pctforeign_controls_placebo" = results_m2_pctforeign_w_controls_placebo)) %>%
        fmtTB(caption="summaries")

rm(placebo_df)

We find the opposite—namely, that the percentage of third-country nationals is, if anything, weakly positively associated with support for far-right parties in previous elections, which casts doubt on the hypothesis that far-right voters are residentially segregated from third-country nationals.

This suggests that the evidence from Vienna specifically for the proposed mechanism—namely, that direct competition with third-country nationals for public housing resources pushes voters to support far-right parties—is weak.

Reweighting

If we had individual rather than aggregate data, we would most likely fit a model at the voter level, regressing whether or not an individual cast a vote for a far-right party in the 2006 federal election against whether they lived in public housing, the number of third-country nationals living in their neighborhood, etc. This paper uses aggregate data at the tract and municipality level as an approximation to the individual data. However, unless we weight the regression by the number of voters in a municipality, we are implicitly weighting voters differently in different municipalities or tracts. As a result, we rerun both models reweighting the observations by the number of voters in the administrative unit.

results_m1_reweight <- estimatr::lm_robust(
  d_rr_06 ~ dv_pop_01 * pct_noneu_06, data = austria_authors.df,
  weights = registered_06)

repro_m1_reweight.df <- tidy_results(list(comp = results_m1_reweight)) %>%
   bind_rows(tidy_m1_original.df)

repro_m1_reweight_summary.df <- tidy_summary(list(comp = results_m1_reweight))

repro_m1_reweight.df %>% group_by(repl_id) %>% 
        fmtTB(caption="coefficients")

repro_m1_reweight_summary.df %>% group_by(repl_id) %>%
      fmtTB(caption="summaries")

We see that both effects become more pronounced—first, that municipalities with higher numbers of third-country residents are substantially less likely to vote for far-right parties, and second, that municipalities with high proportions of people living in public housing are more likely to vote for far-right parties.

We can repeat the same analysis at the level of Vienna tracts for the second analysis.

results_m2_reweight <- 
  estimatr::lm_robust(dv ~ pctrental + pctpublic_w_zsp, data = vienna_authors.df,
                      weights = exp(log_voters), clusters = tract_key)

Warning in eval(quote({: Some observations have missingness in the cluster
variable(s) but not in the outcome or covariates. These observations have been
dropped.

repro_m2_reweight.df <- tidy_results(list(comp = results_m2_reweight)) %>%
  bind_rows(tidy_m2_original.df)

repro_m2_reweight_summary.df <- tidy_summary(list(comp = results_m2_reweight))

repro_m2_reweight.df %>% group_by(repl_id) %>%
        fmtTB(caption="coefficients")

repro_m2_reweight_summary.df %>% group_by(repl_id) %>%  fmtTB(caption="summaries")

In this case, the weighting makes almost no difference to the model results.

We can also perform our alternative specification, where we adjust for the percentage of third-country residents in the Vienna analysis, with reweighting.

results_m2_pctforeign_reweight <-
  estimatr::lm_robust(dv ~ pctforeign * pctpublic_w_zsp, data = vienna_authors.df,
                      weights = exp(log_voters), clusters = tract_key)

tidy_results(list(results_m2_pctforeign_reweight)) %>% fmtTB(caption="coefficients")

tidy_summary(list(results_m2_pctforeign_reweight)) %>% fmtTB(caption="summaries")

Again, this doesn’t seem to meaningfully affect the results.

Data Robustness Analysis

Graphical Exploration of Data (intraoccular impact)

To better understand the authors’ main claims, we attempt to visualize the trends they find in the raw data. We begin with the first analysis, which links support for far-right political parties in the 2006 federal election with the proportion of third-country nationals across Austrian municipalities.

requireNamespace("plotly")

Loading required namespace: plotly

## AUSTRIA
# Plot the relationship between % non-EU and change in vote share
suppressWarnings({
  austria_authors.df %>%
    ggplot(aes(x = pct_noneu_06, y = d_rr_06)) +
    geom_point() +
    labs(x = "% non-EU residents in municipality",
         y = "Change in far-right vote share") +
    # NOTE: There are ~200 municipalities with 0 non-EU residents
    scale_x_log10() +
    geom_smooth(method = lm)
} %>%
  plotly::ggplotly())

`geom_smooth()` using formula = 'y ~ x'

# Plot the relationship between % non-EU and vote share
suppressWarnings({
  austria_authors.df %>%
    ggplot(aes(x = pct_noneu_06, y = rr_share_06)) +
    geom_point() +
    labs(x = "% non-EU residents in municipality",
         y = "Far-right vote share") +
    # NOTE: There are ~200 municipalities with 0 non-EU residents
    scale_x_log10() +
    geom_smooth(method = lm)
} %>%
  plotly::ggplotly())

`geom_smooth()` using formula = 'y ~ x'

# Compare vote share in 2002 and 2006
suppressWarnings({
  austria_authors.df %>%
    select(pct_noneu_06, rr_share_06, rr_share_02) %>%
    pivot_longer(
      cols = starts_with("rr_share"),
      names_prefix = "rr_share_",
      names_to = "year",
      values_to = "vote_share"
    ) %>%
    ggplot(aes(x = pct_noneu_06, y = vote_share)) +
    geom_point() +
    labs(x = "% non-EU residents in municipality",
         y = "Far-right vote share") +
    # NOTE: There are ~200 municipalities with 0 non-EU residents
    scale_x_log10() +
    geom_smooth(method = lm) +
    facet_wrap(vars(year))
} %>%
  plotly::ggplotly())

`geom_smooth()` using formula = 'y ~ x'

Visual examination of the data indicates that (1) there is a very modest association between the percentage of third-country nationals and both the level of support and the change in the level of support for far-right parties in the 2006 elections; and (2) that this trend is more pronounced in 2006 than in 2002. One of this article’s primary interests is, in addition, the interaction between these two factors—i.e., that support is driven by competition between Austrian and third-country nationals for housing, which is most accute when public housing rates and the proportion of the population that is third-country nationals are both high. To visualize this, we plot the relationship between the percentage of third-country nationals and support for far-right parties, stratifying by the percentage of the population that lives in public housing.

# Bin the percentage of the municipality in public housing and plot the
# change in vote share by % non-EU residents
cuts <-
  with(austria_authors.df,
       quantile(dv_pop_01, probs = seq(0, 1, 1 / 3),
                na.rm = TRUE))
suppressWarnings({
  austria_authors.df %>%
    mutate(pct_public_housing = cut(dv_pop_01, breaks = cuts, include.lowest = TRUE)) %>%
    # ~ 13 municipalities don't have data on public housing
    drop_na(pct_public_housing) %>%
    ggplot(aes(x = pct_noneu_06, y = d_rr_06)) +
    geom_point() +
    labs(x = "% non-EU residents in municipality",
         y = "Change in far-right vote share") +
    # NOTE: There are ~200 municipalities with 0 non-EU residents
    scale_x_log10() +
    geom_smooth(method = lm) +
    facet_wrap(vars(pct_public_housing))
} %>%
  plotly::ggplotly())

`geom_smooth()` using formula = 'y ~ x'

We see that while the trend line does appear to get steeper, it is estimated with a considerable amount of uncertainty; in particular, in all three cases, the trend is not visually apparent, and the slope of the fitted linear model has a high degree of uncertainty about sign.

We additionally color the observations based on the baseline level of support for far-right parties to see if an interaction with public housing rates and proportion of third-country nationals is apparent. In particular, it seems plausible that these relationships might be strengthened in places where the baseline support for far-right parties (as measured by their support in 2002) is already high.

# Add coloring based on how far-right the municipality was in the previous
# election
suppressWarnings({
  austria_authors.df %>%
    mutate(pct_public_housing = cut(dv_pop_01, breaks = cuts, include.lowest = TRUE)) %>%
    # ~ 13 municipalities don't have data on public housing
    drop_na(pct_public_housing) %>%
    ggplot(aes(x = pct_noneu_06, y = d_rr_06, color = rr_share_02)) +
    geom_point() +
    labs(x = "% non-EU residents in municipality",
         y = "Change in far-right vote share",
         color = "% far-right in\nprevious election") +
    # NOTE: There are ~200 municipalities with 0 non-EU residents
    scale_x_log10() +
    scale_color_gradient(low = "blue", high = "red") +
    geom_smooth(method = lm) +
    facet_wrap(vars(pct_public_housing))
} %>%
  plotly::ggplotly())

`geom_smooth()` using formula = 'y ~ x'

While some extreme outliers seem to be in places that were already far-right in 2002, no trend is immediately apparent.

Next, we consider the second analysis, which tries to more directly test whether direct competition for public housing between Austrian and third-country nationals explains increasing support for far-right parties in Vienna.

We begin by visualizing the relationship between the proportion of third-country nationals living in a tract, the proportion of residents living in public housing in a tract, and change in far-right vote share.

## VIENNA
# Plot the relationship between % non-EU and change in vote share
suppressWarnings({
  vienna_authors.df %>%
    ggplot(aes(x = pctforeign, y = dv)) +
    geom_point() +
    labs(x = "% non-EU residents in tract",
         y = "Change in far-right vote share") +
    # NOTE: There are ~15 tracts with no foreign residents
    scale_x_log10() +
    geom_smooth(method = lm)
} %>%
  plotly::ggplotly() 
)

`geom_smooth()` using formula = 'y ~ x'

# Plot the relationship between % in public housing and change in vote share
suppressWarnings({
  vienna_authors.df %>%
    ggplot(aes(x = pctpublic_w_zsp, y = dv)) +
    geom_point() +
    labs(x = "% of residents in public housing",
         y = "Change in far-right vote share") +
    # NOTE: There are ~900 tracts with no one in public housing
    scale_x_log10() +
    geom_smooth(method = lm)
  } %>%
  plotly::ggplotly())

`geom_smooth()` using formula = 'y ~ x'

We notice two important facts. First, the relationship between the percentage of individuals living in public housing and the change in far-right vote share is much more pronounced and positive in Vienna than in the national data. Second, the relationship between the percentage of third-country nationals and the far-right vote share is actually very pronounced and negative. To understand the interrelationship between these two factors, we stratify by the proportion of individuals living in public housing and plot the relationship between the percentage of third-country nationals and change in far-right vote share, and see that while the rate of support does seem to rise in the highest bin (i.e., those tracts where the largest number of people live in public housing), the trend within each bin remains fairly negative.

# Plot the relationship between % non-EU and change in vote share, stratified by
# the rate of public housing
cuts <-
  with(vienna_authors.df,
       quantile(
         pctpublic_w_zsp,
         probs = c(0, 1 / 2, 3 / 4, 1),
         na.rm = TRUE
       ))

suppressWarnings({
  vienna_authors.df %>%
    mutate(public_housing = cut(pctpublic_w_zsp, cuts, include.lowest = TRUE)) %>%
    ggplot(aes(x = pctforeign, y = dv)) +
    geom_point() +
    labs(x = "% non-EU residents in tract",
         y = "Change in far-right vote share") +
    # NOTE: There are ~15 tracts with no foreign residents
    scale_x_log10() +
    geom_smooth(method = lm) +
    facet_wrap(vars(public_housing))
} %>%
    plotly::ggplotly())

`geom_smooth()` using formula = 'y ~ x'

rm(cuts)

Outlier analysis

First of all, we can have a look at the outliers in the main variables from Table 1 (Austrian sample). The following code visualizes the relationships between the main variables of interest: the percentage of non-EU residents, the percentage of people living in public housing, and the change in far-right vote share. From each of these variables, we drop 1% of the largest and smallest values. And then the plots compare the change in linear dependence between the variables.

library(patchwork)

###### OUTLIERS PLOTS

p1 <- austria_authors.df %>%
  mutate(outlier = ifelse(
    (
      dv_pop_01 > quantile(dv_pop_01, 0.99, na.rm = T) |
        d_rr_06 > quantile(d_rr_06, 0.99, na.rm = T) |
        dv_pop_01 < quantile(dv_pop_01, 0.01, na.rm = T) |
        d_rr_06 < quantile(d_rr_06, 0.01, na.rm = T)
    ),
    T,
    F
  )) %>%
  ggplot(aes(x = dv_pop_01, y = d_rr_06)) +
  geom_point(aes(color = outlier), size = 1) +
  scale_color_manual(values = c('navy', 'red')) +
  ggtitle('') +
  ylab('Δ  2002–6') + xlab('% public housing') +
  geom_smooth(method = "lm",
              se = T,
              color = "red") +
  geom_smooth(
    data = . %>% filter(outlier == F),
    method = "lm",
    se = T,
    color = "blue"
  )  +
  geom_text(
    data = . %>%
      summarise(r2 = summary(lm(d_rr_06 ~ dv_pop_01))$r.squared),
    aes(
      label = paste("R^2 =", round(r2, 3)),
      x = 0.4,
      y = 0.2
    ),
    color = "red",
    hjust = 0
  ) +
  geom_text(
    data = . %>% filter(outlier == 0) %>%
      summarise(r2 = summary(lm(d_rr_06 ~ dv_pop_01))$r.squared),
    aes(
      label = paste("R^2 =", round(r2, 3)),
      x = 0.4,
      y = 0.18
    ),
    color = "blue",
    hjust = 0
  )

p2 <- austria_authors.df %>%
  mutate(outlier = ifelse(
    (
      pct_noneu_06 > quantile(pct_noneu_06, 0.99, na.rm = T) |
        d_rr_06 > quantile(d_rr_06, 0.99, na.rm = T) |
        pct_noneu_06 < quantile(pct_noneu_06, 0.01, na.rm = T) |
        d_rr_06 < quantile(d_rr_06, 0.01, na.rm = T)
    ),
    T,
    F
  )) %>%
  ggplot(aes(x = pct_noneu_06, y = d_rr_06)) +
  geom_point(aes(color = outlier), size = 1) +
  scale_color_manual(values = c('navy', 'red')) +
  ggtitle('') +
  ylab('Δ  2002–6') + xlab('% non-EU') +
  geom_smooth(method = "lm",
              se = T,
              color = "red") +
  geom_smooth(
    data = . %>% filter(outlier == F),
    method = "lm",
    se = T,
    color = "blue"
  )  +
  geom_text(
    data = . %>%
      summarise(r2 = summary(lm(
        d_rr_06 ~ pct_noneu_06
      ))$r.squared),
    aes(
      label = paste("R^2 =", round(r2, 3)),
      x = 0.15,
      y = 0.2
    ),
    color = "red",
    hjust = 0
  ) +
  geom_text(
    data = . %>% filter(outlier == 0) %>%
      summarise(r2 = summary(lm(
        d_rr_06 ~ pct_noneu_06
      ))$r.squared),
    aes(
      label = paste("R^2 =", round(r2, 3)),
      x = 0.15,
      y = 0.18
    ),
    color = "blue",
    hjust = 0
  )

p3 <- austria_authors.df %>%
  mutate(outlier = ifelse(
    (
      dv_pop_01 > quantile(dv_pop_01, 0.99, na.rm = TRUE) |
        d_rr_06 > quantile(d_rr_06, 0.99, na.rm = TRUE) |
        dv_pop_01 < quantile(dv_pop_01, 0.01, na.rm = TRUE) |
        d_rr_06 < quantile(d_rr_06, 0.01, na.rm = TRUE) |
        pct_noneu_06 > quantile(pct_noneu_06, 0.99, na.rm = TRUE) |
        d_rr_06 > quantile(d_rr_06, 0.99, na.rm = TRUE) |
        pct_noneu_06 < quantile(pct_noneu_06, 0.01, na.rm = TRUE) |
        d_rr_06 < quantile(d_rr_06, 0.01, na.rm = TRUE)
    ),
    TRUE,
    FALSE
  )) %>%
  ggplot(aes(x = I(pct_noneu_06 * pct_noneu_06), y = d_rr_06)) +
  geom_point(aes(color = outlier), size = 1) +
  scale_color_manual(values = c('navy', 'red')) +
  ggtitle('') +  # Add your plot title
  geom_smooth(method = "lm",
              se = T,
              color = "red") + # Add smoothed line for the whole sample
  geom_smooth(
    data = . %>% filter(outlier == F),
    method = "lm",
    se = T,
    color = "blue"
  )  +
  geom_text(
    data = . %>%
      summarise(r2 = summary(lm(
        d_rr_06 ~ I(pct_noneu_06 * pct_noneu_06)
      ))$r.squared),
    aes(
      label = paste("R^2 =", round(r2, 3)),
      x = 0.04,
      y = 0.2
    ),
    color = "red",
    hjust = 0
  ) +
  geom_text(
    data = . %>% filter(outlier == 0) %>%
      summarise(r2 = summary(lm(
        d_rr_06 ~ I(pct_noneu_06 * pct_noneu_06)
      ))$r.squared),
    aes(
      label = paste("R^2 =", round(r2, 3)),
      x = 0.04,
      y = 0.18
    ),
    color = "blue",
    hjust = 0
  )

suppressWarnings(print({p1+p2+p3}))

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

rm(p1,p2,p3)

From these plots, we can observe that after removing outliers, the R-squared values decrease significantly. Respectively, it drops from 0.028 to 0.029 (a 32% decrease), from 0.013 to 0.010 (a 23% decrease), and from 0.008 to 0.003 for the key interaction (a 62.5% decrease).

Now, we replicate the baseline model from Table 1 (regression of far-right vote change on public housing, non-EU and residents and their interaction).

library("lfe")

Loading required package: Matrix

Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack

library("gtsummary")

glance.summary.felm <- function(x,...) {
 x[c("rdf","rse","r.squared","N","r2","r2adj")] %>% as_tibble()
}

austria_authors.df_upd <- austria_authors.df %>%
  mutate(
    outlier98 = ifelse(
      (dv_pop_01 > quantile(dv_pop_01, 0.99, na.rm = TRUE) | 
         d_rr_06 > quantile(d_rr_06, 0.99, na.rm = TRUE) |
         dv_pop_01 < quantile(dv_pop_01, 0.01, na.rm = TRUE) | 
         d_rr_06 < quantile(d_rr_06, 0.01, na.rm = TRUE) |
         pct_noneu_06 > quantile(pct_noneu_06, 0.99, na.rm = TRUE) | 
         d_rr_06 > quantile(d_rr_06, 0.99, na.rm = TRUE) |
         pct_noneu_06 < quantile(pct_noneu_06, 0.01, na.rm = TRUE) | 
         d_rr_06 < quantile(d_rr_06, 0.01, na.rm = TRUE)),
      TRUE, FALSE),
      outlier98_up = ifelse(
        (dv_pop_01 > quantile(dv_pop_01, 0.99, na.rm = TRUE) | 
           pct_noneu_06 > quantile(pct_noneu_06, 0.99, na.rm = TRUE) | 
           d_rr_06 > quantile(d_rr_06, 0.99, na.rm = TRUE)), 
        TRUE, FALSE),
        outlier98_down = ifelse(
      (dv_pop_01 < quantile(dv_pop_01, 0.01, na.rm = TRUE) | 
         pct_noneu_06 < quantile(pct_noneu_06, 0.01, na.rm = TRUE)|
         d_rr_06 < quantile(d_rr_06, 0.01, na.rm = TRUE)), 
      TRUE, FALSE)
    )


m1_ols1.felm <- 
  felm(d_rr_06 ~ dv_pop_01*pct_noneu_06|0|0, data =austria_authors.df)
ols1 <- m1_ols1.felm %>%
  tbl_regression(tidy_fun = purrr::partial(tidy_robust, robust = "HC1"))%>%        
  add_significance_stars(hide_p = F, hide_se = F)%>%
  gtsummary::add_glance_table(include = c(nobs, r.squared))

Arguments `vcov` and `vcov_args` have not been specified in `tidy_robust()`.
Specify at least one to obtain robust standard errors.

tidy_robust(): Robust estimation with
`parameters::model_parameters(model = x, ci = 0.95, robust = "HC1")`

m1_ols2.felm <- 
  felm(d_rr_06 ~ dv_pop_01*pct_noneu_06|0|0, data =subset(austria_authors.df_upd, outlier98 == FALSE))
ols2<- m1_ols2.felm %>%
  tbl_regression()%>%        
  add_significance_stars(hide_p = F, hide_se = F)%>%
  add_glance_table(include = c(nobs, r.squared))  

m1_ols3.felm <-  felm(d_rr_06 ~ dv_pop_01*pct_noneu_06|0|0, data =subset(austria_authors.df_upd, outlier98_up == FALSE)) 

ols3 <-m1_ols3.felm %>%
  tbl_regression()%>%        
  add_significance_stars(hide_p = F, hide_se = F)%>%
  add_glance_table(include = c(nobs, r.squared))


m1_ols4.felm <- felm(d_rr_06 ~ dv_pop_01*pct_noneu_06|0|0, data =subset(austria_authors.df_upd, outlier98_down == FALSE))

ols4<- m1_ols4.felm %>%
  tbl_regression()%>%        
  add_significance_stars(hide_p = F, hide_se = F)%>%
  add_glance_table(include = c(nobs, r.squared))


tbl_merge_ex1 <-
  tbl_merge(
    tbls = list(ols1, ols2, ols3, ols4),
    tab_spanner = c("**Original Specification**", "1% Min-Max Dropped", "1% Max Dropped", "1% Min Dropped")
  )

tbl_merge_ex1

Characteristic	Original Specification			1% Min-Max Dropped			1% Max Dropped			1% Min Dropped
	Beta1	SE2	p-value	Beta1	SE2	p-value	Beta1	SE2	p-value	Beta1	SE2	p-value
dv_pop_01	0.02*	0.012	0.045	0.02	0.014	0.090	0.03	0.015	0.070	0.02	0.012	0.069
pct_noneu_06	-0.02	0.030	0.5	0.02	0.032	0.6	0.02	0.033	0.6	-0.02	0.029	0.4
dv_pop_01 * pct_noneu_06	0.67***	0.183	<0.001	0.47	0.280	0.091	0.50	0.293	0.091	0.67***	0.176	<0.001
No. Obs.	2,373			2,284			2,308			2,349
R²	0.036			0.018			0.019			0.036
1 *p<0.05; **p<0.01; ***p<0.001
2 SE = Standard Error

ml.ls <- list(orig=m1_ols4.felm,maxmintrim=m1_ols4.felm,maxtrim=m1_ols4.felm,mintrim=m1_ols4.felm)
tidy_results(ml.ls) %>%
  fmtTB(caption="coefficients")

tidy_summary(ml.ls)%>%
  fmtTB(caption="summary")

rm(ols1,ols2,ols3,ols4,ml.ls,
   m1_ols4.felm,
   m1_ols4.felm,
   m1_ols4.felm,
   m1_ols4.felm)

Warning in rm(ols1, ols2, ols3, ols4, ml.ls, m1_ols4.felm, m1_ols4.felm, :
object 'm1_ols4.felm' not found

Warning in rm(ols1, ols2, ols3, ols4, ml.ls, m1_ols4.felm, m1_ols4.felm, :
object 'm1_ols4.felm' not found

Warning in rm(ols1, ols2, ols3, ols4, ml.ls, m1_ols4.felm, m1_ols4.felm, :
object 'm1_ols4.felm' not found

We can observe that the initial baseline model demonstrates reproducibility in terms of effect sizes and standard errors. However, once we alter the subsamples, we notice a shift in the results. Specifically, when we exclude the top 1% of the highest values for all three key variables, it leads to a change in the estimate for the interaction term. The coefficient drops to 0.47, with a p-value of 0.091, in contrast to the baseline model where the p-value was less than 0.001.

Moving forward, we explore the impact of clustering errors by district (bezirk). In the initial sample, this adjustment does not lead to a significant alteration (the p-value for the interaction term increases to 0.006). However, it’s noteworthy that when we exclude the highest values, the significance of the term lowers with a p-value of 0.2

m1_rols1.felm <- felm(d_rr_06 ~ dv_pop_01*pct_noneu_06|0|0|bezirk, data =austria_authors.df) 
rols1 <-m1_rols1.felm %>%
  tbl_regression(tidy_fun = purrr::partial(tidy_robust, robust = "HC1"))%>%        
  add_significance_stars(hide_p = F, hide_se = F)%>%
  add_glance_table(include = c(nobs, r.squared))

Arguments `vcov` and `vcov_args` have not been specified in `tidy_robust()`.
Specify at least one to obtain robust standard errors.
tidy_robust(): Robust estimation with
`parameters::model_parameters(model = x, ci = 0.95, robust = "HC1")`

m1_rols2.felm <- felm(d_rr_06 ~ dv_pop_01*pct_noneu_06|0|0|bezirk, data =subset(austria_authors.df_upd, outlier98 == FALSE))
rols2 <-m1_rols2.felm %>%
  tbl_regression()%>%        
  add_significance_stars(hide_p = F, hide_se = F)%>%
  add_glance_table(include = c(nobs, r.squared))  

m1_rols3.felm <- felm(d_rr_06 ~ dv_pop_01*pct_noneu_06|0|0|bezirk, data =subset(austria_authors.df_upd, outlier98_up == FALSE)) 
rols3 <- m1_rols3.felm%>%
  tbl_regression()%>%        
  add_significance_stars(hide_p = F, hide_se = F)%>%
  add_glance_table(include = c(nobs, r.squared))

m1_rols4.felm <- felm(d_rr_06 ~ dv_pop_01*pct_noneu_06|0|0|bezirk, data =subset(austria_authors.df_upd, outlier98_down == FALSE)) 
rols4 <- m1_rols2.felm %>% 
    tbl_regression()%>%        
  add_significance_stars(hide_p = F, hide_se = F)%>%
  add_glance_table(include = c(nobs, r.squared))


tbl_merge_ex1 <-
  tbl_merge(
    tbls = list(rols1, rols2, rols3, rols4),
    tab_spanner = c("**Original Specification**", "1% Min-Max Dropped", "1% Max Dropped", "1% Min Dropped")
  )

tbl_merge_ex1

Characteristic	Original Specification			1% Min-Max Dropped			1% Max Dropped			1% Min Dropped
	Beta1	SE2	p-value	Beta1	SE2	p-value	Beta1	SE2	p-value	Beta1	SE2	p-value
dv_pop_01	0.02	0.024	0.3	0.02	0.028	0.4	0.03	0.029	0.4	0.02	0.028	0.4
pct_noneu_06	-0.02	0.060	0.7	0.02	0.064	0.8	0.02	0.065	0.8	0.02	0.064	0.8
dv_pop_01 * pct_noneu_06	0.67**	0.242	0.006	0.47	0.373	0.2	0.50	0.374	0.2	0.47	0.373	0.2
No. Obs.	2,373			2,284			2,308			2,284
R²	0.036			0.018			0.019			0.018
1 *p<0.05; **p<0.01; ***p<0.001
2 SE = Standard Error

ml.ls <- list(orig=m1_rols4.felm, 
              maxmintrim=m1_rols3.felm, 
              maxtrim=m1_rols3.felm,
              mintrim=m1_rols4.felm)

tidy_results(ml.ls) %>%
  fmtTB(caption="coefficients")

tidy_summary(ml.ls)%>%
  fmtTB(caption="summary")

rm(ml.ls, rols1, rols2, rols3, rols4, m1rols1.felm, m1rols2.felm,m1rols3.felm,m1rols4.felm)

Warning in rm(ml.ls, rols1, rols2, rols3, rols4, m1rols1.felm, m1rols2.felm, :
object 'm1rols1.felm' not found

Warning in rm(ml.ls, rols1, rols2, rols3, rols4, m1rols1.felm, m1rols2.felm, :
object 'm1rols2.felm' not found

Warning in rm(ml.ls, rols1, rols2, rols3, rols4, m1rols1.felm, m1rols2.felm, :
object 'm1rols3.felm' not found

Warning in rm(ml.ls, rols1, rols2, rols3, rols4, m1rols1.felm, m1rols2.felm, :
object 'm1rols4.felm' not found

Here we look at the covariate balance between outliers and non-outlier observations.

library("cobalt")

 cobalt (Version 4.5.1, Build Date: 2023-04-27)

library("MatchIt")

Attaching package: 'MatchIt'

The following object is masked from 'package:cobalt':

    lalonde

austria_authors.df_upd1 <- subset(austria_authors.df_upd, outlier98_up == TRUE | outlier98_up == FALSE)



austria_authors.df_upd1 <- austria_authors.df_upd1 %>% filter(!is.na(educ_tertiary) & !is.na(avg_income) &
                                !is.na(lab_pct_manufact_01) & !is.na(lab_pct_unemp) &
                                !is.na(welfare_cap_06 ) & !is.na( health_cap_06) &
                                !is.na(education_cap_06) & !is.na(foreignborn_delta))

set.seed(7)
m.out <- MatchIt::matchit(outlier98_up ~ educ_tertiary + registered_06+
                            avg_income + 
                            lab_pct_manufact_01 + 
                            lab_pct_unemp + welfare_cap_06 +  
                            health_cap_06 + education_cap_06 + 
                            foreignborn_delta, data = austria_authors.df_upd1)

bal.tab(m.out, thresholds = c(m = .1), un = TRUE)

Balance Measures
                        Type Diff.Un Diff.Adj        M.Threshold
distance            Distance  0.6520   0.3392                   
educ_tertiary        Contin.  0.4443   0.1043 Not Balanced, >0.1
registered_06        Contin.  0.5805   0.3800 Not Balanced, >0.1
avg_income           Contin.  0.2420   0.0743     Balanced, <0.1
lab_pct_manufact_01  Contin. -0.0888  -0.1974 Not Balanced, >0.1
lab_pct_unemp        Contin.  0.8857   0.1788 Not Balanced, >0.1
welfare_cap_06       Contin.  0.6576   0.2847 Not Balanced, >0.1
health_cap_06        Contin.  0.5589   0.3480 Not Balanced, >0.1
education_cap_06     Contin.  0.6096   0.2370 Not Balanced, >0.1
foreignborn_delta    Contin.  0.6280   0.2128 Not Balanced, >0.1

Balance tally for mean differences
                   count
Balanced, <0.1         1
Not Balanced, >0.1     8

Variable with the greatest mean difference
      Variable Diff.Adj        M.Threshold
 registered_06     0.38 Not Balanced, >0.1

Sample sizes
          Control Treated
All          2304      65
Matched        65      65
Unmatched    2239       0

love.plot(m.out, stats = c("mean.diffs"),
          thresholds = c(m = .1, v = 2), abs = TRUE, 
          binary = "std",
          var.order = "unadjusted")

rm(m.out,austria_authors.df_upd1)

Finally, we apply the baseline specification to the sample of 65 outliers. In this case, we observe a significantly higher R-squared value (0.687). Moreover, all the coefficients are highly significant with p-values less than 0.001, and they exhibit comparatively large effect sizes. For instance, the coefficient of the key interaction is 3.0, which is 4.5 times higher than that of the overall sample. This suggests that the observed effect in the full sample in the paper may be largely driven by this small subsample.

We can observe that the sample of outlier observations (depicted in red on the plot) significantly differs from the rest in several aspects. Parameters such as the number of registered voters, education, income, and various other factors tend to be higher on average in this subgroup. It’s worth noting that this smaller sample also exhibits higher rates of unemployment, increased welfare and education spending, as well as a more pronounced growth in the number of foreign-born residents.

out1 <- felm(d_rr_06 ~ dv_pop_01*pct_noneu_06|0|0|bezirk, data =subset(austria_authors.df_upd, outlier98_up == TRUE)) %>%
  tbl_regression()%>%        
  add_significance_stars(hide_p = F, hide_se = F)%>%
  add_glance_table(include = c(nobs, r.squared))

tbl_merge_ex2 <-
  tbl_merge(
    tbls = list(out1),
    tab_spanner = c("")
  )
tbl_merge_ex2 

Characteristic
	Beta1	SE2	p-value
dv_pop_01	-0.26***	0.035	<0.001
pct_noneu_06	-0.98***	0.069	<0.001
dv_pop_01 * pct_noneu_06	3.0***	0.332	<0.001
No. Obs.	65
R²	0.687
1 *p<0.05; **p<0.01; ***p<0.001
2 SE = Standard Error

Model fit

sqe1 <- double(nrow(austria_authors.df))
for (i in seq(nrow(austria_authors.df))) {
  row <- austria_authors.df[i, ]
  m <- lm(d_rr_06 ~ dv_pop_01 * pct_noneu_06, data = austria_authors.df[-i, ])
  sqe1[[i]] <- (row$d_rr_06 - predict(m, row))^2
}
print(glue::glue("LOO RMSE for interaction mode: { sqrt(mean(sqe1, na.rm = TRUE)) }"))

LOO RMSE for interaction mode: 0.0283408162629289

sqe2 <- double(nrow(austria_authors.df))
for (i in seq(nrow(austria_authors.df))) {
  row <- austria_authors.df[i, ]
  m <- lm(d_rr_06 ~ dv_pop_01 + pct_noneu_06, data = austria_authors.df[-i, ])
  sqe2[[i]] <- (row$d_rr_06 - predict(m, row))^2
}
print(glue::glue("LOO RMSE for main effects: { sqrt(mean(sqe2, na.rm = TRUE)) }"))

LOO RMSE for main effects: 0.0284132384077463

rm(sqe1)
rm(sqe2)

sqe1 <- double(nrow(vienna_authors.df))
for (i in seq(nrow(vienna_authors.df))) {
  row <- vienna_authors.df[i, ]
  m <- lm(dv ~ pctrental * pctpublic_w_zsp, data = vienna_authors.df[-i, ])
  sqe1[[i]] <- (row$dv - predict(m, row))^2
}
print(glue::glue("LOO RMSE for interaction mode: { sqrt(mean(sqe1, na.rm = TRUE)) }"))

LOO RMSE for interaction mode: 0.0431238824470698

sqe2 <- double(nrow(vienna_authors.df))
for (i in seq(nrow(vienna_authors.df))) {
  row <- vienna_authors.df[i, ]
  m <- lm(dv ~ pctrental + pctpublic_w_zsp, data = vienna_authors.df[-i, ])
  sqe2[[i]] <- (row$dv - predict(m, row))^2
}
print(glue::glue("LOO RMSE for main effects: { sqrt(mean(sqe2, na.rm = TRUE)) }"))

LOO RMSE for main effects: 0.0439284411572522

rm(sqe1)
rm(sqe2)

As an addenda, after reflection upon the authors’ responses to the first version of this comment, we provide further clarification on how the outlier analysis should be interpreted. The results analysis provides evidence that the data, as a whole does not conform to the favored model specification ( a linear model with a simple interaction term ). It is not intended as a definitive rebuttal of the conceptual theory advanced in the authors’ publication – but to highlight the need for more extensive attention to model diagnostics and to the need for explicit comparison with competing models.

Conclusion

We find that the research is computationally replicable but the evidence for the main causal claims are overstated. While the focus of analysis concentrates on patterns in outcomes that are likely relevant to understanding the underlying data-generating process, the statistical models supporting the causal claims support only a small proportion of overall variance in outcomes. Furthermore, the analysis elides relevant competing models.

The relative weakness of the claim is obscured in the public analysis because neither overall goodness of fit, nor comparison to ‘naive’ / baseline models are included. While conceptual reproducibility analysis would be useful for exploring more complex alternate models this route is obstructed by the absence of citations, documentation and linking codes that would support reliable reanalysis using original data, or augmentation of the authors’ data with additional measures. We conjecture that, as a general practice, research reliability would be increased by including these practices in publication and data sharing.

References

Altman, Micah, Jeff Gill, and Michael McDonald. 2004. Numerical issues in statistical computing for the social scientist. Wiley series in probability and statistics. Hoboken, NJ: Wiley-Interscience.

Cavaille, Charlotte. 2023. “Replication Data for: “How Distributional Conflict over in-Kind Benefits Generates Support for Far-Right Parties.” Harvard Dataverse. https://doi.org/10.7910/DVN/SYNP73.

Cavaillé, Charlotte, and Jeremy Ferwerda. 2023. “How Distributional Conflict over In-Kind Benefits Generates Support for Far-Right Parties.” The Journal of Politics 85 (1): 19–33. https://doi.org/10.1086/720643.

Kadane, Joseph B, and Nicole A Lazar. 2004. “Methods and Criteria for Model Selection.” Journal of the American Statistical Association 99 (465): 279–90. https://doi.org/10.1198/016214504000000269.

Footnotes

1. Corresponding author.

2. All authors declare that they have no financial support or conflict of interest in this publication.

3. See https://osf.io/3dseb